// so here's an idea... only perturb one parameter, at random...


struct Rn *perturb_eigen_Rn_1 (struct Rn *data) { // THIS CAN BE USER-SPECIFIED
        struct Rn *return_data=cp_Rn(data);

	return_data->f[(uint)u01()*(data->n)]+=fabs(STDEV*n01(u01(),u01(),0));//fabs(gsl_ran_gaussian(rg,STDEV));

        return return_data;
}


struct Rn *perturb_Rn_1 (struct Rn *data,struct svd *ed) { // THIS CAN BE USER-SPECIFIED
	
	struct Rmn *covariance=theorize_covariance_Rmn(ed);

	uint rand=u01()*data->n;

	struct Rn *return_data=cp_Rn(data);

	return_data->f[rand]+=n01(u01(),u01(),0)*covariance->f[rand]->f[rand];

        return return_data;
}

struct Rmn *perturb_Rmn_1 (struct Rmn *data) {
	struct Rmn *return_data=cp_Rmn(data);	 //copy our initial data...

	return_data->f[(uint)u01()*(data->m)]->f[(uint)u01()*(data->n)]+=STDEV*n01(u01(),u01(),0);//gsl_ran_gaussian(rg,STDEV);

	return return_data;
}

struct Rmn *perturb_Rmn (struct Rmn *data) {
	struct Rmn *return_data=init_Rmn(data->m,data->n,1);	 //generate a random PSD covariance matrix

	for (uint i=0;i<data->m;i++)
		for (uint j=0;j<data->n;j++)
			return_data->f[i]->f[j]+=data->f[i]->f[j]; //add to our existing one.

	return return_data;
}
struct Rn *sample_mean (struct mpi_progress d) { // this only works on mean...
	return d.mean->f[(uint)u01()*(mpi_ranks*accepts_per_rank)];
}

float new_threshold(struct mpi_progress d) {

//	mergeSort(d.rejection_threshold,ACCEPTANCES);
	float f=0;
	for (uint i=0;i<ACCEPTANCES;i++)
		f+=d.rejection_threshold[i];
	f/=(1.f*ACCEPTANCES); //just taknig the average threshold...

//	cerr<<"low="<<d.rejection_threshold[0]<<",median="<<d.rejection_threshold[ACCEPTANCES/2]<<",mean="<<f<<",high="<<d.rejection_threshold[ACCEPTANCES-1]<<endl;


	return f;//d.rejection_threshold[ACCEPTANCES/2];
}


struct Rn *perturb_Rn (struct Rn *data, struct svd *ed, uint df) { // THIS CAN BE USER-SPECIFIED

	struct Rmn *perturbed_Rmn=theorize_Rmn(data, ed, df);

	struct Rn *return_data=mean_Rn(perturbed_Rmn);

	destroy_Rmn(perturbed_Rmn);

	//float u1,u2;
	//for (uint i=0;i<data->n;i++) {
	//	if (i%2==0)
	//		u1=u01(), u2=u01();
	//	return_data->f[i]+=STDEV*n01(u1,u2,i%2);//gsl_ran_gaussian(rg,STDEV);
//	}
	return return_data;
}

/*
struct Rmn *theorize_Rmn (struct Rn *mean, struct Rmn *covariance,uint m) { // THIS IS THE HEART OF IT....
	assert(mean->n==covariance->n);
	uint n=mean->n;
	struct Rmn *return_data=init_Rmn(m,n,0);
	double *d1=NULL,*d2=NULL,*d3=NULL;

	for (uint i=0;i<m;i++) {
		destroy_Rn(return_data->f[i]);
		return_data->f[i]=array_to_Rn(d1=ranMVNormal(d2=Rmn_to_array(covariance),d3=Rn_to_array(mean),n),n);
		free(d1), free(d2), free(d3);
	}

	return return_data;

}*/
struct Rmn_summary *summarize_Rmn (struct Rmn *a){
	struct Rmn_summary *d=init_Rmn_summary(a->n);

	d->x1=moment_Rn(a,1);
	d->x2=moment_Rn(a,2);


	return d;
}


float sumsquared_Rn(struct Rn *a, struct Rn *b) {
assert(a->n==b->n);

float sum=0;
for (uint i=0;i<a->n;i++)
	sum+=powf(a->f[i]-b->f[i],2);

return sum;


}


float sumsquared_Rmn(struct Rmn *a, struct Rmn *b) {
assert(a->m==b->m && a->n==b->n);

float sum=0;
for (uint i=0;i<a->m;i++)
	sum+=sumsquared_Rn(a->f[i],b->f[i]);
	//for (uint j=0;j<a->n;j++)
		//sum+=pow(a->f[i]->f[j]-b->f[i]->f[j],2);
return sum;
}


float exchangeable_sumsquared_Rmn(struct Rmn *a, struct Rmn *b) {

assert(a->m==b->m && a->n==b->n);


float eps,tmp,sum=0;
char disqualified[a->m];
for (uint i=0;i<a->m;i++)
	disqualified[i]=0;
uint eps_n=0;
for (uint i=0;i<a->m;i++) {
	eps=FLT_MAX;

	for (uint j=0;j<a->m;j++)
		if (!disqualified[j]) //if not disqualified
			if ((tmp=sumsquared_Rn(a->f[i],b->f[j]))<eps)
				eps_n=j, eps=tmp;
	//disqualify:
	disqualified[eps_n]=1;
	//add to sum
	sum+=eps;
}

return sum;



}

struct Rmn *normalize_Rmn (struct Rmn *d) {  // we are using MoM here, which may need a sampling bias adjustment on the variance...


struct Rmn *r=init_Rmn(d->m,d->n,0);

struct Rn *x1=moment_Rn(d,1);
struct Rn *x2=moment_Rn(d,2);


for (uint i=0;i<d->m;i++)
	for (uint j=0;j<d->n;j++) {
		r->f[i]->f[j]=(d->f[i]->f[j]-x1->f[j])/(sqrt(x2->f[j]-x1->f[j]*x1->f[j]));

	}


return r;

}

double kolmogorov_p(double x,uint n) { //gives TWO SAMPLE p-value from critical value...

double r=0,d;
x/=sqrt(2.f*n/(1.f*n*n));


for (uint k=1;(d=exp(-1.f*(2.f*k-1)*(2.f*k-1)*PI*PI/(8.f*x*x)))>0;k++) 
	r+=d;

return 1-sqrt(2.f*PI)/x*r;

}

double kolmogorov_smirnov_compare(struct Rn *d1, struct Rn *d2) { //returns a p-value that these come from the same distribution

assert(d1->n == d2->n);
mergeSort(d1->f,d1->n);
mergeSort(d2->f,d2->n);

double x=0, sup_x=0;

	for (uint i=0;i<d1->n;i++) { //testing d1->f[i] as our X_i... tryintg to find how many under this
		for (uint j=0;j<d1->n;j++) { //should this be j=i ?
			if (d1->f[i] < d2->f[j]) {
				x=j/(double)d2->n;
				break;
			}

		}
	if (fabs((i+1)/(double)d1->n-x)>sup_x)
		sup_x=fabs((i+1)/(double)d1->n-x);
	}

return sup_x;

}

double distance_Rmn(struct Rmn *a, struct Rmn *b) {
/*	float f=0;

#define MOMENTS	2 //must be at least two in order to standardize
	assert(a->m==b->m && a->n==b->n);

	//struct Rmn_summary *a_summary=summarize_Rmn(a), *b_summary=summarize_Rmn(b);
	struct Rmn *a_summary=init_Rmn(MOMENTS,a->n,0), *b_summary=init_Rmn(MOMENTS,b->n,0);

	for (uint i=0;i<MOMENTS;i++) {
		destroy_Rn(a_summary->f[i]);
		destroy_Rn(b_summary->f[i]);
		a_summary->f[i]=moment_Rn(a,i+1);
		b_summary->f[i]=moment_Rn(b,i+1);
	}

	for (uint i=0;i<MOMENTS;i++) {		
		for (uint j=0;j<a->n;j++) {
//			cerr<<i<<"moment is: "<<a_summary->f[i]->f[j]<<"and standardized moment is:"<<a_summary->f[i]->f[j]/powf(sqrt(a_summary->f[1]->f[j]-a_summary->f[0]->f[j]*a_summary->f[0]->f[j]),i+1)<<endl;

			f+=powf(a_summary->f[i]->f[j]/powf(sqrt(a_summary->f[1]->f[j]-a_summary->f[0]->f[j]*a_summary->f[0]->f[j]),i+1)-b_summary->f[i]->f[j]/powf(sqrt(b_summary->f[1]->f[j]-b_summary->f[0]->f[j]*b_summary->f[0]->f[j]),i+1),2);

//			f+=powf(a_summary->f[i]->f[j]/powf((a_summary->f[1]->f[j]-a_summary->f[0]->f[j]*a_summary->f[0]->f[j]),i+1)-b_summary->f[i]->f[j]/powf((b_summary->f[1]->f[j]-b_summary->f[0]->f[j]*b_summary->f[0]->f[j]),i+1),2); // this is supposedly a standardized moment...
		//	f+=powf(a_summary->f[i]->f[j]-b_summary->f[i]->f[j],2);
		}
	}

	destroy_Rmn(a_summary);
	destroy_Rmn(b_summary);


	return sqrt(f);
*/
//	return sumsquared_Rmn(a,b);//exchangeable_sumsquared_Rmn(a,b);
	//write_Rmn(a,cerr);
	//cerr<<"vs"<<endl;
//	write_Rmn(b,cerr);
	
	//transpose...


if (METRIC==TSS)
	return sumsquared_Rmn(a,b);
else if (METRIC==KS){
	struct Rn *a_n=init_Rn(a->m,0), *b_n=init_Rn(a->m,0);
	for (uint i=0;i<a->m;i++) {
		a_n->f[i]=a->f[i]->f[0];
		b_n->f[i]=b->f[i]->f[0];
	}
double x=kolmogorov_p(kolmogorov_smirnov_compare(a_n,b_n),a->m);
	destroy_Rn(a_n), destroy_Rn(b_n);
	return 1-x;
}

}
